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This is the second in a series of papers on the construction and validation of a three-dimensional code for the 
solution of the coupled system of the Einstein equations and of the general relativistic hydrodynamic equations, 
and on the application of this code to problems in general relativistic astrophysics. In particular, we report on the 
accuracy of our code in the long-term dynamical evolution of relativistic stars and on some new physics results 
obtained in the process of code testing. The following aspects of our code have been validated: the generation 
of initial data representing perturbed general relativistic polytropic models (both rotating and non-rotating), the 
long-term evolution of relativistic stellar models, and the coupling of our evolution code to analysis modules 
providing, for instance, the detection of apparent horizons or the extraction of gravitational waveforms. The 
tests involve single non-rotating stars in stable equilibrium, non-rotating stars undergoing radial and quadrupo- 
lar oscillations, non-rotating stars on the unstable branch of the equilibrium configurations migrating to the 
stable branch, non-rotating stars undergoing gravitational collapse to a black hole, and rapidly rotating stars 
in stable equilibrium and undergoing quasi-radial oscillations. We have carried out evolutions in full general 
relativity and compared the results to those obtained either with perturbation techniques, or with lower dimen- 
sional numerical codes, or in the Cowling approximation (in which all the perturbations of the spacetime are 
neglected). In all cases an excellent agreement has been found. The numerical evolutions have been carried out 
using different types of polytropic equations of state using either the rest-mass density only, or the rest-mass 
density and the internal energy as independent variables. New variants of the spacetime evolution and new high 
resolution shock capturing (HRSC) treatments based on Riemann solvers and slope limiters have been imple- 
mented and the results compared with those obtained from previous methods. In particular, we have found the 
"monotonized central differencing" (MC) limiter to be particularly effective in evolving the relativistic stellar 
models considered. Finally, we have obtained the first eigenfrequencies of rotating stars in full general relativity 
and rapid rotation. A long standing problem, such frequencies have not been obtained by other methods. Over- 
all, and to the best of our knowledge, the results presented in this paper represent the most accurate long-term 
three-dimensional evolutions of relativistic stars available to date. 



I. INTRODUCTION 

Computational general relativistic astrophysics is an in- 
creasingly important field of research. Its development is be- 
ing driven by a number of factors. Firstly, the large amount of 
observational data by high-energy X-ray and 7-ray satellites 
such as Chandra, XMM and others Secondly, the new 
generation of gravitational wave detectors coming online in 
the next few years [^], and thirdly, the rapid increase in com- 
puting power through massively parallel supercomputers and 
the associated advance in software technologies, which make 
large-scale, multi-dimensional numerical simulations possi- 
ble. Three-dimensional (3D) simulations of general relativis- 
tic astrophysical events such as stellar gravitational collapse 
or collisions of compact stars and black holes are needed to 



fully understand the incoming wealth of observations from 
high-energy astronomy and gravitational wave astronomy. It 
is thus not surprising that in recent years hydrodynamical sim- 
ulations of compact objects in numerical relativity has become 
the focus of several research groups ||3 - 1 1 . 

In a previous paper [|[) (hereafter paper I) we presented 
a 3D general-relativistic hydrodynamics code (GR_Astro) 
constructed for the NASA Neutron Star Grand Challenge 
Project [pH. The GR_Astro code has been developed by 
Washington University and the Albert Einstein Institute and 
has the capability of solving the coupled set of the Einstein 
equations and the general relativistic hydrodynamic (GRHy- 
dro) equations [p^]. It has been built using the Cactus Com- 
putational Toolkit constructed by the Albert Einstein In- 
stitute, Washington University and other institutes. Paper I 
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presented our formulation for the GRHydro equations coupled 
either to the standard Arnowitt-Deser-Misner (ADM) [0] for- 
mulation of the Einstein equations or to a hyperbolic formu- 
lation of the equations [|l5|]. It demonstrated the consistency 
and convergence of the code for a comprehensive sample of 
testbeds having analytic solutions. It gave a detailed analy- 
sis of twelve different combinations of spacetime and hydro- 
dynamics evolution methods, including Roe's and other ap- 
proximate Riemann solvers, as well as their relative perfor- 
mance and comparisons when applied to the various testbeds. 
The code as described and validated in paper I has been ap- 
plied to various physical problems, such as those discussed in 
Refs. [^[jT|,|l^, and is now freely available [|T|]. 

The main purpose of this paper is to examine and validate 
our code in long-term, accurate simulations of the dynamics 
of isolated stars in strong gravitational fields. Single relativis- 
tic stars are indeed expected as the end-point of a number of 
astrophysical scenarios (such as gravitational collapse and bi- 
nary neutron star merging) and should provide important in- 
formation about strong field physics both through electromag- 
netic and gravitational wave emissions. A number of new nu- 
merical techniques have been incorporated in the present code 
leading to a much improved ability to simulate relativistic 
stars. These techniques concern both the evolution of the field 
equations, for which we have implemented new conformal- 
traceless formulations of the Einstein equations, and in the 
evolution of the hydrodynamical variables, for which the use 
of the "monotonized central differencing" (MC) limiter has 
provided us with the small error growth-rates necessary for 
simulations over several dynamical timescales. 

More precisely, in this paper we focus on the accuracy of 
the code during long-term evolution of spherical and rapidly 
rotating stellar models. We also investigate the nonlinear dy- 
namics of stellar models that are unstable to the fundamen- 
tal radial mode of pulsation. Upon perturbation, the unstable 
models will either collapse to a black hole, or migrate to a con- 
figuration in the stable branch of equilibrium configurations 
(a behavior studied in the case of unstable boson stars ijl^]). 
In the case of collapse, we follow the evolution all the way 
down to the formation of a black hole, tracking the genera- 
tion of its apparent horizon. In the case of migration to the 
stable branch, on the other hand, we are able to accurately 
follow the nonlinear oscillations that accompany this process 
and that can give rise to strong shocks. The ability to simu- 
late large amplitude oscillations is important as we expect a 
neutron star formed in a supernova core-collapse [ p^|pO| ] or 
in the accretion-induced collapse of a white dwarf to oscillate 
violently in its early stages of life. 

Particularly important for their astrophysical implications, 
we study the Unear pulsations of spherical and rapidly rotating 
stars. The computed frequencies of radial, quasi-radial and 
quadrupolar oscillations are compared with the corresponding 
frequencies obtained with lower-dimensional numerical codes 
or with alternative techniques such as the Cowling approx- 
imation (in which the spacetime is held fixed and only the 
GRHydro equations are evolved) or relativistic perturbative 
methods. The comparison shows an excellent agreement con- 
firming the ability of the code to extract physically relevant 



information from tiny perturbations. The successful determi- 
nation of the eigenfrequencies for rapidly rotating stars com- 
puted with our code is noteworthy. Such frequencies have not 
been obtained before with the system being too complicated 
for perturbative techniques. 

The simulations discussed here make use of two different 
polytropic equations of state (EOS). In addition to the stan- 
dard "adiabatic" EOS, in which the pressure is expressed as 
a power law of the rest-mass density, we have carried out 
simulations implementing the "ideal fluid" EOS, in which the 
pressure is proportional to both the rest-mass density and the 
specific internal energy density. This latter choice increases 
the computational costs (there is one additional equation to 
be solved) but allows for the modeling of non-adiabatic pro- 
cesses, such as strong shocks and the conversion of bulk ki- 
netic energy into internal energy, which are expected to ac- 
company relativistic astrophysical events. 

There are a number of reasons why we advocate the careful 
validation of general relativistic astrophysics codes. Firstly, 
the space of solutions of the coupled system of the Ein- 
stein and GRHydro equations is, to a large extent, unknown. 
Secondly, the numerical codes must solve a complicated set 
of coupled partial differential equations involving thousands 
of terms and there are plenty of chances for coding errors. 
Thirdly, the complex computational infrastructure needed for 
the use of the code in a massively parallel environment, in- 
creases the risk of computational errors, a risk that can only be 
minimized through meticulous tests such as those presented 
here as well as in paper I. This paper, however, wants to be 
more than a list of testbeds: the results presented show that 
our current numerical methods are mature enough for obtain- 
ing answers to new and outstanding problems in the physics 
of relativistic stars. 

The organization of this paper is as follows: the formulation 
of the differential equations for the spacetime and the hydro- 
dynamics is briefly reviewed in section H Section ^ gives a 
short description of the numerical methods, with emphases on 
the new schemes introduced in this paper (in addition to those 
in paper I). Sections |iv| - |v| represent the core of the paper 
and there we present the main results of our simulations. In 
section |^ in particular, we focus our attention on the simula- 
tion of nonrotating relativistic stars. In section^ we consider 
the evolution of rotating stars. Section VI is dedicated to the 
extraction of gravitational waveforms generated by the non- 
radial pulsations of perturbed relativistic stars. In section VII 
we summarize our results and conclusions. We use a spacelike 
signature (— , +, +, +) and units in which c = G — Mq — 1 
(geometric units based on solar mass) unless explicitly spec- 
ified. Greek indices are taken to run from to 3 and Latin 
indices from 1 to 3. 



II. BASIC EQUATIONS 

We give a brief overview of the system of equations in this 
section. We refer the reader to paper I for more details. 
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A. Field Equations 

In general relativity, the dynamics of the spacetime is de- 
scribed by the Einstein field equations G^^ — SttT^^, with 
Gfj.1, being the Einstein tensor and Tf^^, the stress-energy ten- 
sor. Many different formulations of the equations have been 
proposed throughout the years, starting with the ADM for- 
mulation in 1962 [[l4|]. In our code, we have implemented 
three different formulations of the field equations, includ- 
ing the ADM formulation, a hyperbolic formulation [|l5| and 
a more recent conformal-traceless formulation based on the 
ADM construction [|l]j22|| (see also Ref. [|3ll). 

In the ADM formulation [p^], the spacetime is foliated with 
a set of non-intersecting spacelike hypersurfaces. Two kine- 
matic variables relate the surfaces: the lapse function a, which 
describes the rate of advance of time along a timelike unit vec- 
tor n'^ normal to a surface, and the shift three-vector /3* that 
relates the spatial coordinates of two surfaces. In this con- 
struction the line element reads 

ds^ = -{a^ - (3iP')df + 2(3idx'dt + j.jdx^dx^ . (1) 

The original ADM formulation casts the Einstein equations 
into a first-order (in time) quasi-linear [^] system of equa- 
tions. The dependent variables are the 3-metric jij and the 
extrinsic curvature Kij . The evolution equations read 



-8^ (Su 



(2) 



Rj 



(3) 



where V; denotes the covariant derivative with respect to the 
3-metric 7^, i?y is the Ricci curvature of the 3-metric, and 
K = j'^-' Kij is the trace of the extrinsic curvature. In addition 
to the evolution equations, there are four constraint equations: 
the Hamiltonian constraint 



and the momentum constraints 



(4) 



Vj K'^ - 7'^ K - Stt/ = 0. (5) 

In equations (@)-(|l), Padm^/' ^tj^^ = 1^^ are the com- 
ponents of the stress-energy tensor projected onto the 3D hy- 
persurface (for a more detailed discussion, see Ref. [^). 

As mentioned above, in addition to the two formulations de- 
scribed in paper I, we have recently implemented a conformal- 
traceless reformulation of the ADM system, as proposed 
by [ pTp^ ]. Details of our particular implementation of this 
formulation are extensively described in Ref. | p3| ] and will 
not be repeated here. We only mention here that this for- 
mulation makes use of a conformal decomposition of the 3- 



metric, 7^ 
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■fij and the trace-free part of the extrin- 



sic curvature, Aij — Kij — ^ijK/i, with the conformal fac- 
tor (j) chosen to satisfy e"*"^ — 7^/'^ = det(7y )^/'^. In this 



formulation, as shown in Ref. [^, in addition to the evo- 
lution equations for the conformal three-metric 7^^ and the 
conformal-traceless extrinsic curvature variables Aij, there 
are evolution equations for the conformal factor (f>, the trace 
of the extrinsic curvature K and the "conformal connection 
functions" f ' (following the notation of Ref. [p2[]). We note 
that the final mixed, first and second-order, evolution system 
for 10, K, 'jij, Aij,r^} is not in any immediate sense hyper- 
bolic [P6|]. In the original formulation of Ref. [|T]|, the auxil- 
iary variables Fi = — J^j were used instead of the r\ 



In Refs. [23 27[| the improved properties of this conformal- 
traceless formulation of the Einstein equations were compared 
to the ADM system. In particular, in Ref. a number of 
strongly gravitating systems were analyzed numerically with 
convergent HRSC methods with total-variation-diminishing 
(TVD) schemes using the equations described in paper I. 
These included weak and strong gravitational waves, black 
holes, boson stars and relativistic stars. The results show 
that our treatment leads to a stable numerical evolution of the 
many strongly gravitating systems. However, we have also 
found that the conformal-traceless formulation requires grid 
resolutions higher than the ones needed in the ADM formula- 
tion to achieve the same accuracy, when the foliation is made 
using the "K-driver" approach discussed in Ref. [^]. Because 
in long-term evolutions a small error growth-rate is the most 
desirable property, we have adopted the conformal-traceless 
formulation as our standard form for the evolution of the field 
equations. 



B. Hydrodynamic Equations 

The GRHydro equations are obtained from the local conser- 
vation laws of the density current (continuity equation) and of 
the stress-energy tensor, which we assume to be that of a per- 
fect fluid T^"' = phut^u" + Pgt"", with u^' being the fluid 
4-velocity and P and h the (isotropic) pressure and the spe- 
cific enthalpy, respectively. In our code the GRHydro equa- 
tions are written as a first-order flux-conservative hyperbolic 
system [|9||] 



dtU + d^F^ = S , 



(6) 



where the evolved state vector U is given, in terms of the 
primitive variables: the rest-mass density p, the 3-velocity 
— jW + fi"^ ja and the specific internal energy e, as 



(7) 



Here 7 is the determinant of the 3-metric 7^^ and W is the 
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Lorentz factor, W — avP = (1 — jijv'^v^) 
more, the 3-flux vectors F^ are given by 



-1/2 



Further- 
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F''= a[{v'^^P')S, + ^P5^^ . (8) 

Finally, the source vector S is given by 



5= a^T>^-g,„T%, , (9) 

where F"^^ are the Christoffel symbols. 

C. Gauge Conditions 

The code is designed to handle arbitrary shift and lapse con- 
ditions, which can be chosen as appropriate for a given space- 
time simulation. More information about the possible families 
of spacetime slicings which have been tested and used with 
the present code can be found in Refs. [^^. Here, we limit 
ourselves to recall details about the specific foliations used in 
the present evolutions. In particular, we have used algebraic 
slicing conditions of the form 



{dt - [fdi)a = -f{a) a^K, 



(10) 



with f{a) > but otherwise arbitrary. This choice contains 
many well known slicing conditions. For example, setting 
/ = 1 we recover the "harmonic" slicing condition, or by 
setting f = q/a, with q being an integer, we recover the gen- 
eralized "iH-log" slicing condition [ ^0[ ] which after integration 
becomes 



a = g{x') + -log7, 



(11) 



where g(a;*) is an arbitrary function of space only. In partic- 
ular, all of the simulations discussed in this paper are done 
using condition ( pi] ) with q = 2, basically due to its compu- 
tational efficiency (we caution that "gauge pathologies" could 
develop with the "l-nlog" slicings, see Ref. [ pTp2[ |). 

The evolutions presented in this paper were carried out with 
the shift vector being either zero or constant in time. 



III. NUMERICAL METHODS 

We now briefly describe the numerical schemes used in our 
code. We will distinguish the schemes implemented in the 
evolution of the Einstein equations from those implemented 
in the evolution of the hydrodynamic equations. In both cases, 
the equations are finite-differenced on spacelike hypersurfaces 
covered with 3D numerical grids using Cartesian coordinates. 



A. Spacetime Evolution 

As described in paper I, our code supports the use of sev- 
eral different numerical schemes [||,^. Currently, a Leapfrog 
(non-staggered in time) and an iterative Crank-Nicholson 
scheme have been coupled to the hydrodynamic solver. 

The Leapfrog method assumes that all variables exist on 
both the current time step t" and the previous time step 
Variables are updated from to (future time) evalu- 
ating all terms in the evolution equations on the current time 
step t". The iterative Crank-Nicholson solver, on the other 
hand, first evolves the data from the current time step to 
the future time step using a forward in time, centered in 
space first-order method. The solution at steps t" and 
are then averaged to obtain the solution on the half time step 
^n+i/2 -pjjjg solution at the half time step is then used 

in a Leapfrog step to re-update the solution at the final time 
step This process is then iterated. The error is defined as 
the difference between the current and previous solutions on 
the half time step This error is summed over all grid- 

points and all evolved variables. Because the smallest number 
of iterations for which the iterative Crank-Nicholson evolution 
scheme is stable is three and further iterations do not improve 
the order of convergence [ ^3l^3| ], we do not iterate more than 
three times. Unless otherwise noted, all simulations reported 
in this paper use the iterative Crank-Nicholson scheme for the 
time evolution of the spacetime. 



B. Hydrodynamical Evolution 

The numerical integration of the GRHydro equations 
is based on High-Resolution Shock-Capturing (HRSC) 
schemes, specifically designed to solve nonlinear hyperbolic 
systems of conservation laws. These conservative schemes 
rely on the characteristic structure of the equations in order to 
build approximate Riemann solvers. In paper I we presented 
a spectral decomposition of the GRHydro equations, suitable 
for a general spacetime metric (see also Ref. [Ql). 

Approximate Riemann solvers compute, at every cell- 
interface of the numerical grid, the solution of local Riemann 
problems (i.e. the simplest initial value problem with dis- 
continuous initial data). Hence HRSC schemes automatically 
guarantee that physical discontinuities developing in the solu- 
tion (e.g., shock waves, which appear in core-collapse super- 
novae or in coalescing neutron star binaries) are treated con- 
sistently. HRSC schemes surpass traditional approaches [|^,||] 
which rely on the use of artificial viscosity to resolve such dis- 
continuities, especially for large Lorentz factor flows. HRSC 
schemes have a high order of accuracy, typically second-order 
or more, except at shocks and extremal points. We refer the 
reader to [ ^ , pq ] for recent reviews on the use HRSC schemes 
in relativistic hydrodynamics. 

One of the major advantages of using HRSC schemes is 
that we can take advantage of the many different algorithms 
that have been developed and tested in Newtonian hydrody- 
namics. In this spirit, our code allows for three alternative 
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ways of performing the numerical integration of the hydrody- 
namic equations: (i) using a flux-split method [^; (ii) using 
Roe's approximate Riemann solver [^], and (Hi) using Mar- 
quina's flux-formula [|3^]. The different methods differ simply 
in the way the numerical fluxes at the cell-interfaces are calcu- 
lated in the corresponding flux-formula. The code uses slope- 
limiter methods to construct second-order TVD schemes 
by means of monotonic piecewise linear reconstructions of the 
cell-centered quantities to the left (L) and right (R) sides of ev- 
ery cell-interface for the computation of the numerical fluxes. 
More precisely, and are computed to second-order 
accuracy as follows: 



- Ui + diix^^i - Xi) 



Xi+l) 



(12) 
(13) 



where x denotes a generic spatial coordinate. We have fo- 
cused our attention on two different types of slope limiters, 
the standard "minmod" limiter and the "monotonized central- 
difference" (MC) limiter [pi]]. In the first case, the slope Ci is 
computed according to 



fjj — minmod 



Ax 



Ax 



(14) 



where Aa; denotes the cell spacing. The minmod function of 
two arguments is defined by 

a if \a\ < \b\ and ab > 
minmod(a,5) = if |fe| < \a\ and ab > 
if a6 < 

On the other hand, the MC slope limiter (which was not in- 
cluded in the previous version of the code discussed in paper 
I) does not reduce the slope as severely as minmod near a dis- 
continuity and, therefore, a sharper resolution can be obtained. 
In this case the slope is computed as 



MC 



Ax 



Ax 



(15) 



where the MC function of two arguments is defined by 

' 2a if \a\ < and 2\a\ < |c|, and ab > 

2b if |6| < |a|, and 2\b\ < |c|, and ab > 
c if |c| < 2|a|, and |c| < 2|5|, anda5 > 
if a6 < 



MC(a,6) 



and where c = {a + b) /2. Both schemes provide the desired 
second-order accuracy for smooth solutions, while still satis- 
fying the TVD property. In sect. IV A we will report on a 
comparison between the two algorithms and justify the use of 
the MC slope limiter as our preferred one. 



C. Equations of State 

As mentioned in the Introduction, to explore the behavior 
of our code in long-term evolutions of equilibrium configura- 
tions, we used two different polytropic equations of state and 
at various central rest-mass densities. In particular, we have 
implemented both an adiabatic (or zero temperature) EOS 



l+l/JV 



P ^Kp^ = Kp 
and as a so-called "ideal fluid" EOS 

P = (r - l)pe , 



(16) 



(17) 



where K is the polytropic constant, F the polytropic index 
and N = {V ~ 1)^^ the polytropic exponent. The ideal fluid 
EOS (|l7|) depends on both the rest-mass density p and on the 
specific internal energy e, it corresponds to allowing the poly- 
tropic coefficient K in adiabatic EOS ( p^ to be a function of 
entropy. The use of an adiabatic EOS with a constant K is 
computationally less expensive and is physically reasonable 
when modeling configurations that are in near equilibrium, 
such as stable stellar models in quasi-equilibrium evolutions. 
There are however dynamical processes, such as those involv- 
ing nonlinear oscillations and shocks, in which the variations 
in the energy entro py can not be neglected. The simulations 
discussed in section [V C , where both equations of state (|l6|)- 
( [l^ ) are used for the same configuration, gives direct evidence 
of how a more realistic treatment of the internal energy of the 
system can produce qualitatively different results. 

The increased accuracy in the physical description of the 
dynamical system comes with a non-negligible additional 
computational cost. It involves the solution of an additional 
equation (i.e. the evolution equation for the specific internal 
energy e), increasing the total number of GRHydro equations 
from four to five and making accurate long-term evolutions 
considerably harder. 



D. Boundary Conditions 

In our general-purpose code, a number of different bound- 
ary conditions can be imposed for either the spacetime vari- 
ables or for the hydrodynamical variables. We refer the reader 
to ||^,p3[| for details. In all of the runs presented in this paper 
we have used static boundary conditions for the hydrodynam- 
ical variables and radiative outgoing boundary conditions for 
the spacetime variables. The only exception to this is the evo- 
lution of rotating stars (see sect. for which the spacetime 
variables have also been held fixed at the outer boundary. 



IV. SPHERICAL RELATIVISTIC STARS 

We turn next to the description of the numerical evolutions 
of relativistic star configurations. We start by considering 
spherical models. 
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A. Long-term evolution of stable configurations 

Using isotropic coordinates {t, r, 9, </>), the metric describ- 
ing a static, spherically symmetric relativistic star reads 



2 ■ 2 

r sm 



'), (18) 



where v and A are functions of the radial coordinate r only. 
The form of the metric component g^r is much simpler in 
these coordinates than in Schwarzschild coordinates, which 
are often used to describe a Tolman-Oppenheimer-VoUcoff 
(TOV) equilibrium stellar solution. In addition, is not con- 
strained to be equal to unity at the center of the stellar config- 
uration, as in Schwarzschild coordinates. We have found that 
these two properties of the isotropic coordinates are very ben- 
eficial to achieve long-term numerical evolutions of relativis- 
tic stars. Therefore, all simulations of spherical relativistic 
stars shown in this paper have been performed adopting the 
line element (|l8|) expressed in Cartesian coordinates. 
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FIG. 1. Evolution of the central rest-mass density pc (in units of 
the initial central rest-mass density Pc,o) for a nonrotating star with 
gravitational mass M = 1.65 Mq. Using Roe's approximate Rie- 
mann solver, the figure shows a comparison in the use of the minmod 
and of the MC slope limiters for both the ideal fluid and the adiabatic 
EOS. 

Although the initial configurations refer to stellar models 
in stable equilibrium, the truncation errors at the center and 
at the surface of the star excite small radial pulsations that 
are damped in time by the numerical viscosity of the code. 
Moreover, these pulsations are accompanied by a secular evo- 
lution of the values of the central rest-mass density away 
from its initial value. Similar features have been reported in 
Refs. | p^ , p3[ |. These features converge away at the correct rate 
with increasing grid resolution and do not influence the long- 
term evolutions. Moreover, the secular evolution of the central 
rest-mass density varies according to the EOS adopted: when 



using the ideal fluid EOS, we have observed that the secular 
drift of the central rest-mass density is towards lower densi- 
ties. However, if we enforced the adiabatic condition (which 
is justified for the case of a near-equilibrium evolution), we 
have observed that the dominant truncation error has opposite 
sign and the central rest-mass density evolves towards larger 
values. 



2.00 



a 




1.00 



0.75 

3 

t(ms) 

FIG. 2. Evolution of the normalized central rest-mass density pc 
for a nonrotating M — 1.65 Mq star. Different lines show a com- 
parison between Roe's Riemann solver and Marquina's flux-formula 
for different slope limiters. 

This is shown in Fig. |l] where we plot the evolution of a 
TOV star with gravitational mass M — 1.65 Mq, constructed 
with a iV = 1 polytrope. In our units, the polytropic con- 
stant is K = 123.5 and the central rest-mass density of the 
star is pc = 1.00 x 10^'^. For these tests, a very coarse grid 
of 35"^ gridpoints in octant symmetry is sufficient and allows 
the major effects to be revealed with minimal computational 
costs. The outer boundary is placed at about 1.7 (where r^. 
is the isotropic coordinate radius of the star). We use radiative 
boundary conditions with a l/r fall-off. Irrespective of the 
slope limiter used, the magnitude of the secular drift observed 
in the central rest-mass density evolution is roughly a factor 
of two smaller for the adiabatic EOS than for the ideal fluid 
EOS. As a result, in all of the evolutions of stable configu- 
rations which remain close to equilibrium (such as pulsating 
stars, with no shock developing), the adiabatic EOS is pre- 
ferred. 

Fig. [l| also gives a comparison of the use of the minmod 
and the MC slope limiters in the evolution of the normalized 
central rest-mass density. For both the ideal fluid and the adia- 
batic EOS, the MC limiter shows a significantly smaller secu- 
lar increase in the central rest-mass density, as compared to the 
minmod one. The simulations in Fig. |l] employed Roe's ap- 
proximate Riemann solver in the fluid evolution scheme and 
this is then compared to Marquina's flux-formula in Fig. § 
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for the evolution of the central rest-mass density. The secular 
increase is significantly smaller when using Marquina's flux- 
formula than when using Roe's solver, and this is especially 
noticeable for the minmod slope limiter. A comparison of the 
increase of the maximum error in the Hamiltonian constraint 
after several ms of evolution (not shown here) indicates that it 
is about 80% smaller with Marquina than with Roe, when us- 
ing the adiabatic EOS. As a result of the above comparisons, 
we have adopted Marquina's scheme with the MC slope lim- 
iter as our preferred scheme for evolution of the GRHydro 
equations. Unless otherwise noted, all of the simulations pre- 
sented in this paper have been obtained with such a scheme. 
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FIG. 3. Time evolution of the normalized central rest-mass den- 
sity at three different grid resolutions (32"', 64'^ and 96"^ gridpoints, 
respectively), for a M — 1.4 Mq, N = 1 relativistic, spherical 
polytrope. The evolution of the central rest-mass density is mainly 
modulated by the fundamental radial mode of oscillation of the star. 
The initial amplitude of the oscillation converges to zero at sec- 
ond-order, while the secular increase in the central rest-mass density 
converges away at almost second-order. 

Next, we show in Fig. ^ the long-term evolution of the 
central rest-mass density for three different grid resolutions. 
For this, we consider a nonrotating N = 1 polytropic star 
with gravitational mass M = 1.4 Mq, circumferential radius 
R — 14.15 km, central rest-mass density pc — 1.28 x 10"'^ 
and K — 100. The different simulations used 32'^, 64'^ 
and 96"^ gridpoints with octant symmetry and with the outer 
boundary placed at 1.7 Vg. These grid resolutions correspond 
to about 19, 38 and 56 gridpoints per star radius, respectively. 
Fig. ^ shows the oscillations in the central rest-mass density 
and the secular evolution away from the initial value men- 
tioned above. The oscillations are produced by the first-order 
truncation error at the center and the surface of the star (our 
hydrodynamical evolution schemes are globally second order, 
but only first-order at local extrema; see related discussions in 
Ref. [P3[|, where long-term convergence tests are presented) 



but both the amplitude of the initial oscillation and the rate 
of the secular change converge to zero at nearly second-order 
with increasing grid resolution. 

Note that the evolutions shown in Figs. |^-|| extend to 7 ms, 
corresponding to about 10 dynamical times (taking the fun- 
damental radial mode period of pulsation as a measure of the 
dynamical timescale), significantly longer than, say, the ones 
reported by other authors [|8}0| . Our evolutions are limited 
by the time available (a simulation with 96^ gridpoints and up 
to 7 ms takes about 40 hours on a 128 processor Cray-T3E su- 
percomputer.). We have found that for a resolution of 96'^, the 
central density at the end of the 7ms evolution is just 0.25% 
larger than the initial central density. 

For the same configuration, we show, in Fig. ^ the time 
evolution of the L2-norm of the violation of the Hamiltonian 
constraint at the three different grid resolutions. Also in this 
case, the violation of the Hamiltonian constraint converges to 
zero at nearly second-order with increasing grid resolution. 
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FIG. 4. Convergence of the L2-norm of the Hamiltonian con- 
straint, at three different grid resolutions (32'^, 64'' and 96'^ grid- 
points, respectively), for a i\f = 1.4 Mq, N — 1 polytropic spheri- 
cal relativistic star The rate of convergence is close to second-order 
with increasing grid resolution. 

In Fig. I], we show other aspects of the accuracy of the sim- 
ulation with 96^ gridpoints, by comparing the initial profiles 
of the rest-mass density p and of the lapse function a of the 
TOV star with those obtained after 7 ms of evolution. The 
small deviations from the original profiles are worth empha- 
sizing. The small inset shows a magnification of the rapid 
change in the gradient of the rest-mass density profile at the 
surface of the star 
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FIG. 5. Variation of the original profiles along the x-axis of the 
rest-mass density (left vertical axis) and lapse function (right vertical 
axis), for a M =1.4 Mq, N = 1 polytropic spherical relativistic 
star, after 7 ms of evolution. A 96^ grid in octant symmetry was used 
in the simulation. The small inset shows a magnification of the rapid 
change in the gradient of the rest-mass density profile at the surface 
of the star. 



B. Radial pulsations 

As mentioned in the previous section, the truncation errors 
of the hydrodynamical schemes used in our code trigger radial 
pulsations of the initially static relativistic star (see Ref. [ p3| ] 
for a review). These pulsations are initiated at the surface of 
the star, where the gradients of the rest-mass density are the 
largest (cf. Fig. [s]). Because gravitational waves cannot be 
emitted through the excitation of radial pulsations of nonrotat- 
ing relativistic stars, these pulsations are damped only by the 
numerical viscosity of the code in numerical simulations of 
inviscid stars. In treatments more dissipative than the HRSC 
schemes used in our code, such as those using artificial vis- 
cosity or particle methods (e.g. Smoothed Particle Hydrody- 
namics), these oscillations will be damped significantly faster. 

In order to test the properties of the long-term hydrodynam- 
ical evolution separately from those of the spacetime evolu- 
tion, we have first examined the long-term hydrodynamical 
evolution separately from those of the spacetime evolution, we 
have first examined the small-amplitude radial pulsations in a 
fixed spacetime of an initially static relativistic star. As initial 
data, we use the M = 1.4 Mq polytropic star of the previous 
section. We show, in Fig. ||, the evolution up to 7 ms of the 
normalized star's central rest-mass density with a numerical 
grid of 96'^ gridpoints. The amplitude of the excited pulsa- 
tions in this purely hydrodynamical evolution is minute (less 
than 1 part in 200) and is significantly smaller than the corre- 
sponding amplitude in a coupled hydrodynamical and space- 
time evolution (compare the vertical axes of Figs. I and |). 
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FIG. 6. Time evolution of the central rest-mass density of a 
A/ = 1.4 A/0, A'^ = 1 polytropic spherical relativistic star. In 
this the simulation the spacetime is held fixed and the hydrodynamic 
variables have been evolved on a numerical grid of 96"^ gridpoints. 
The evolution is a superposition of radial normal modes of pulsation, 
excited by truncation errors of the hydrodynamical scheme. Higher 
overtones are damped faster by the small but non-zero numerical vis- 
cosity. 



A closer look at Figure ^ reveals that the evolution of the 
central rest-mass density is a superposition of different radial 
normal modes of pulsation. The higher-frequency modes are 
damped faster, so that after a certain time the evolution pro- 
ceeds mainly in the fundamental mode of pulsation. Note also 
the small damping rate of the fundamental pulsation mode, in- 
dicating the small effective numerical viscosity of our HRSC 
hydrodynamical scheme. The evolution towards larger values 
of the central rest-mass density is similar to that discussed in 
Section IV A but less pronounced in this case. At a resolution 
of 96"^ gridpoints, the secular change in the average central 
rest-mass density is less than 0.02% for the total evolution 
time shown. 

The use of truncation error as an initial perturbation de- 
serves commenting on. The oscillations caused by truncation 
error will converge away with increasing resolution, hence the 
overall oscillation amplitude can carry no physical informa- 
tion about the system. However, the frequencies and normal- 
ized eigenfuntions of particular normal-modes of oscillation 
of the star are physical (in the sense that they match the eigen- 
frequencies and eigenfunctions calculated through perturba- 
tive analyses) and can be extracted from these simulations 
by carrying out a Fourier transform of the time evolution of 
the radial velocity or of the rest-mass density. As the small- 
amplitude pulsations are in the linear regime, the eigenfunc- 
tions can be normalized arbitrarily (e.g. to 1 .0 at the surface of 
the star). At increasing resolution, the solution converges to 
the mode-frequencies and to the normalized eigenfunctions. 
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even though the overall oscillation amplitude converges to 
zero. Such evolutions are useful for extracting the properties 
of linear normal-modes of oscillation, as long as the resolu- 
tion is fine enough that the pulsations excited by truncation 
errors are in the linear regime and as long as the resolution is 
coarse enough that the various local 1st and 2nd order trunca- 
tion errors of the numerical scheme result in a time evolution 
that is dominated by a sum of normal modes (at very fine res- 
olutions the Fourier transform of the time evolution would be 
very small and thus have a very noisy power spectrum due to 
roundoff errors, in which case the physical normal-mode fre- 
quencies would be difficult to extract - this has not been the 
case for the resolutions used in this paper). We also note that 
different variants of our hydrodynamical evolution schemes 
excite the various physical normal-modes at different ampli- 
tudes. For example, 2nd order schemes employing the min- 
mod limiter tend to clearly excite a large number of high- 
frequency overtones, whereas the use of the MC limiter results 
in the clear excitation of only a few low-frequency overtones 
and a more noisy FFT power spectrum at higher frequencies 
(for the resolutions used in this paper). This difference in be- 
haviour is due to the differences in the local truncation errors 
inherent in these numerical schemes. 

The radial pulsations are a sum of eigen modes of pulsa- 
tion. Since the radial pulsations triggered by truncation errors 
have a small amplitude, one can compare the frequencies with 
that computed by linear perturbation theory or with hy- 
drodynamical evolutions of similar models in 2D [ |4^ , p3| ]. In 
this way we can validate that the "artificial" perturbations pro- 
duced by the truncation errors do excite "physical" modes of 
oscillation for a relativistic star However, before discussing 
the results of this comparison, it is important to emphasize 
that the identification of the frequency peaks in the Fourier 
transform of the time evolution of a given variable with phys- 
ical frequencies must be done with care. A real pulsation fre- 
quency must be global (the same at every point in the star, at 
least for discrete normal mode frequencies) and it should ap- 
pear in the time evolution of different physical quantities de- 
scribing the star's structure and dynamics. To eliminate possi- 
ble ambiguities, we have carried out our frequency identifica- 
tion procedure for different variables and at different positions 
in the star 

Fig. shows the Fourier transform of the time evolution of 
the central rest-mass density of the same initial model as in 
Fig. I, but using the minmod limiter (which gives a clearer ex- 
citation of the higher overtones). We indicate with F the fun- 
damental normal mode frequency and with HI — H6 the next 
six higher frequency modes (overtones). We have also com- 
pared the frequency peaks in the Fourier spectrum to both the 
normal mode frequencies expected by linear perturbation the- 
ory in the Cowling approximation (see Ref. [p6|]), and to the 
frequencies computed with an independent 2D axisymmetric 
nonlinear code [^3]], which uses the same HRSC schemes but 
in spherical polar coordinates (shown as dashed vertical lines 
in Fig. 0). 

As can be seen from Table |, the agreement is extremely 
good. The relative difference between the 3D and 2D results 
at this grid resolution is better than 1% up to (HA) and slightly 




f (kHz) 

FIG. 7. Fourier transform of the central rest-mass density evolu- 
tion of a Af = 1.4 Mq, N — 1 polytropic spherical relativistic star, 
in a fixed spacetime evolution. Here F represents the fundamental 
normal mode frequency, while HI — H6 indicate the first six over- 
tones. The frequency peaks in the power spectrum are in excellent 
agreement with the radial normal mode frequencies (shown here as 
dashed vertical lines) computed with an independent 2D code using 
spherical polar coordinates. The solid and dotted lines were com- 
puted with 96'^ and 64^ gridpoints, respectively. The units of the 
vertical axis are arbitrary. 



larger for higher frequencies (H5 and H6), which become 
under-resolved at this grid resolution. This excellent agree- 
ment is a significant test for the correct implementation of the 
hydrodynamical evolution schemes in our code, and is an indi- 
cation of the level of accuracy we can achieve, resolving and 
following these small deviations away from the equilibrium 
configuration. As one would expect, lower or higher resolu- 
tion runs (e.g. with 64'^ or 144'^ gridpoints), which have intrin- 
sically larger or smaller perturbation amplitudes, respectively, 
reproduce the peaks in the power spectrum shown in Fig. ^ 
(see dotted line in Fig. ^ which corresponds to an evolution 
with 64'^ grid-points. 

After establishing the accuracy of the long-term evolution 
of the GRHydro equations, we have examined the eigenfre- 
quencies of the radial pulsations of spherical stars in coupled 
hydrodynamical and spacetime evolutions. A Fourier trans- 
form of the evolution of the radial velocity (for the same equi- 
librium model as the one discussed before) is shown in Fig. |[ 
Again in this case, we have been able to identify several fre- 
quency peaks in the Fourier spectrum with the normal mode 
frequencies obtained with linear perturbation techniques [|47|]. 
A detailed comparison of these frequencies is shown in Ta- 
ble p. The agreement is again excellent. Note the rather large 
differences between the frequencies shown in Tables | and |^ 
The Cowling approximation is rather inaccurate for the lowest 
radial mode-frequencies [Esl], but is increasingly more accu- 
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FIG. 8. Fourier transform of the evolution of the radial velocity 
for a M = 1.4 A'Iq, N = 1 polytropic spherical relativistic star in 
a coupled spacetime and hydrodynamical evolution. The frequency 
peaks in the spectrum are in excellent agreement with the radial nor- 
mal mode frequencies computed by perturbation theory (shown here 
as dashed vertical lines). As in Fig. M, here F represents the funda- 
mental normal mode frequency, while HI — H3 are the next three 
higher frequency modes. The units of the vertical axis are arbitrary. 

rate for nonradial pulsations or for higher frequencies [p^]. 

All of the results discussed so far refer to simulations in- 
volving stable relativistic configurations. In the following 
section we consider numerical evolutions of relativistic stars 
which are initially in an unstable equilibrium. 

C. Migration of unstable configurations to the stable branch 

The numerical evolution of a nonrotating, relativistic star 
in an equilibrium unstable to the fundamental radial mode of 
pulsation is mainly determined by the numerical truncation 
errors that cause it to evolve away from its initial configura- 
tion. Depending on the type of perturbation, the star can either 
collapse to a black hole or expand and migrate to the stable 
branch of the sequence of equilibrium models, reaching a new, 
stable equilibrium configuration with approximately the same 
rest-mass of the perturbed star We have therefore constructed 
a model of a TV = 1, A' = 100 polytropic star with rest-mass 
Mo = 1.535 Mq (M = 1.447 Af©) and a central rest-mass 
density pc = 8.0 x 10~^, which is larger than the central rest- 
mass density of the maximum-mass stable model. The star is 
therefore initially in an unstable equilibrium (see the inset of 
Fig. ^ and under the perturbation introduced by the truncation 
error, it expands, evolving rapidly to smaller central rest-mass 
densities, until it reaches the stable branch of equilibrium con- 
figurations. An analogous behavior has been observed in nu- 
merical simulations of relativistic boson stars [ p^ (see also 



TABLE 1. Comparison of small-amplitude radial pulsation fre- 
quencies obtained with the present 3D nonlinear evolution code with 
frequencies obtained with an independent 2D code. Both codes 
evolve the GRHydro equations in a fixed spacetime and for an equi- 
librium model of a A*' = 1 relativistic polytrope with AI/R = 0.15. 



Mode 


Present 3D code 
(kHz) 


2D code 
(kHz) 


Relative Difference 

(%) 


F 


2.696 


2.701 


0.2 


HI 


4.534 


4.563 


0.6 


H2 


6.346 


6.352 


0.1 


H3 


8.161 


8.129 


0.4 


H4 


9.971 


9.875 


1.0 


H5 


11.806 


11.657 


1.3 


H6 


13.605 


13.421 


1.7 



TABLE II. Comparison of small-amplitude radial pulsation fre- 
quencies obtained with the present 3D nonlinear evolution code with 
linear perturbation mode frequencies, in fully coupled evolutions. 
The equilibrium model is a nonrotating A'^ = 1 relativistic polytrope 
with A//i? = 0.15. 



Mode 


Present 3D code 
(kHz) 


Perturbation code 
(kHz) 


Relative Difference 

(%) 


F 


1.450 


1.442 


0.6 


HI 


3.958 


3.955 


0.0 


H2 


5.935 


5.916 


0.3 


H3 


7.812 


7.776 


0.4 



Ref. for recent numerical simulations of expanding un- 
stable boson stars). 

In a realistic astrophysical scenario, a stable neutron star 
can accrete matter e.g. from a companion star in a binary sys- 
tem or from infalling matter after its formation in a supernova 
core-collapse. The star would then secularly move towards 
larger central densities along the stable branch of equilibrium 
configurations, exceed the maximum-mass limit and collapse 
to a black hole. No secular mechanism could evolve the star 
to the unstable branch. In this respect, the migration mech- 
anism discussed here cannot occur in practice. Nevertheless, 
it provides a consistent solution of the initial value problem 
and represents an important test of the accuracy of the code 
in a highly dynamical and non-adiabatic evolution. We use 
such an initial data set to study large amplitude oscillations 
of relativistic stars, which cannot be treated accurately by lin- 
ear perturbation theory. Large amplitude oscillations about a 
configuration on the stable branch could occur after a super- 
nova core-collapse [ pO[ ] or after an accretion-induced collapse 
of a white dwarf. While the actual set of quasi-normal modes 
excited will depend on the excitation process, the ability to 
simulate large amplitude oscillations is important. 

Fig. ^ shows the evolution of the central rest-mass density 
Pc normalized to its initial value and up to a final time of 4.26 
ms. On a very short dynamical timescale of 0.5 ms the star 
has expanded and has its central density dropped to about 3 
% of its initial central rest-mass density. Note that this is less 
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FIG. 9. Evolution of the (normalized) central rest-mass density pc 
during the migration of an unstable relativistic star to a stable model 
with the same rest-mass. When an adiabatic EOS is used (dotted 
line) the difference in gravitational binding energy between the un- 
stable and stable models is periodically converted in bulk kinetic en- 
ergy through highly nonlinear, nearly constant amplitude pulsations. 
In contrast, when an ideal fluid EOS is used (solid line), the grav- 
itational binding energy is gradually converted into internal energy 
via shock heating. As a result, the oscillations are damped and the 
heated stable equilibrium model approaches a central density slightly 
smaller than the rest-mass density of a zero temperature star of the 
same rest-mass (indicated by an asterisk on the left vertical axis). 



than the central rest-mass density, pc = 1-35 x 10~^, of the 
stable model of same rest-mass, which is indicated with an as- 
terisk on the vertical axis of Fig. ^ During the rapid decrease 
of the central rest-mass density, the star acquires a large radial 
momentum. The star then enters a phase of large amplitude 
radial oscillations about the stable equilibrium model with the 
same rest-mass. Because the unstable and stable models have 
rather different degrees of compactness, the migration to the 
stable branch will be accompanied by the release of a signif- 
icant amount of gravitational binding energy which could ei- 
ther be converted to bulk kinetic energy or to internal energy 
depending on the choice of EOS. 

In order to investigate both responses, we have performed 
two different evolutions of the same initial model. In the first 
case (the "adiabatic EOS" in Fig. H), we have enforced the adi- 
abatic condition during the evolution, i.e. we have assumed 
that the star remains at zero temperature following an adia- 
batic EOS. As shown in Fig. ^ with a dotted line, in this case 
the star behaves like a compressed spring which is allowed 
to expand, oscillating with a nearly constant amplitude. This 
indicates that the star periodically converts all of the excess 
gravitational binding energy into the kinetic energy and vice 
versa. As the oscillations are highly nonlinear, the restoring 
force is weaker at higher densities than at lower densities and 
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FIG. 10. Shock formation in the outer core mantle, during the first 
bounce at equilibrium densities of an unstable star, evolved with an 
ideal fluid EOS. The top and bottom panels show the internal energy 
e and radial velocity v,^, respectively, at three different times: the 
homologous infall phase, the inner core bounce and the outwards 
shock propagation. The oscillations of the inner core are damped by 
shock heating. 



the oscillations are therefore far from being sinusoidal. 

In the second case (the "ideal fluid EOS" in Fig. we do 
not enforce the abiabatic condition, but allow all of thermody- 
namic variables to evolve in time. As a result, the oscillations 
are gradually damped in time, while the star oscillates around 
a central density close to that of a stable star with the same 
rest-mass. 

The rapid decrease in the oscillation amplitude is due to the 
dissipation of kinetic energy via shock heating. At the end of 
the first expansion (i.e. at the first minimum in Fig. the 
star has expanded almost to the edge of the numerical grid. At 
this point, the outer parts of the initial star have formed a low- 
density, outer-core mantle around the high-density inner core 
and the star then starts to contract. Fig. |l^ shows with solid 
lines the supersonic infall of the outer core mantle at < = 0.84 
ms, while the inner core is contracting homologously. Af- 
ter this "point of last good homology", the high-density inner 
core reaches its maximum infall velocity and then starts slow- 
ing down. The infalling low-density mantle forms a shock at 
the inner core/mantle boundary (dotted lines at t = 0.98 ms in 
Fig. |l^. After the inner core bounces, it expands and pressure 
waves at the inner core-mantle boundary feed the shock wave 
with kinetic energy (dashed lines at i = 1.13 ms in Fig. [l0| ). 
In this way, the shock wave is dissipating the initial binding 
energy of the star so that the amplitude of the central density 
oscillations decreases with time. The above process is very 
similar to the core bounce in neutron star formation (see, for 
instance, the description in [50|), except for the fact that here 
the outer mantle is created during the first rapid expansion 
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from material of the initial unstable star. 

As a result of the damping of the radial oscillations, the star 
settles down, on a secular timescale, to a stable equilibrium 
configuration with central density somewhat smaller than the 
central density of a stable star with same rest-mass as the ini- 
tial unstable star. This is because part of the matter of the 
initial star forms a heated mantle around the inner core. 

The evolution shown in Fig. ^ was obtained using a resolu- 
tion of 96^ gridpoints. Since the initial unstable configuration 
is much more compact than the final configuration, the bound- 
aries of the computational grid were placed at about 4.5 r^. As 
a result, the grid resolution of the initial configuration is rather 
low, causing an additional, non-negligible deviation of the av- 
erage central rest-mass density of the pulsating star away from 
the expected central rest-mass density of the zero-temperature 
star of the same rest-mass. 

The evolution of the highly nonlinear and nonadiabatic pul- 
sations of a star when it settles down on the stable branch, un- 
derlines the importance of evolving all of the thermodynamic 
variables (including the specific internal energy) and the im- 
portance of using HRSC methods in order to resolve the for- 
mation and evolution of shocks correctly. These capabilities 
of the numerical code will be important in the correct simu- 
lation of general relativistic astrophysical events such as the 
merging of a neutron star binary system or the formation of a 
neutron star in an accretion-induced collapse of a white dwarf. 



D. Gravitational collapse of unstable configurations 

As mentioned in the previous section, the numerical scheme 
used in the hydrodynamical evolution is such that it causes a 
nonrotating relativistic star in an unstable equilibrium to ex- 
pand and migrate to the configuration of same rest-mass lo- 
cated on the stable branch of equilibrium configurations. In 
order to study the gravitational collapse to a black hole of an 
unstable model we need to add to the initial model a small ra- 
dial perturbation in the rest-mass density distribution. A very 
small perturbation of the order of '--^ 1% is sufficient and its ra- 
dial dependence can be simply given by cos(7rr/2rs), where 
r is coordinate distance from the center and its value at 
the surface of the star. The addition of this small perturba- 
tion dominates over the truncation error and causes the star 
to collapse to a black hole. Note that after the perturbation is 
added to the initial equilibrium configuration, the constraint 
equations are solved to provide initial data which is a solution 
to the field equations . 

The (forced) collapse to a black hole of an unstable spher- 
ical relativistic star is shown in Fig. [Ill for a simulation with 
128'^ gridpoints in octant symmetry, using Roe's solver and an 
ideal fluid EOS. The figure shows the profiles along the x-axis 
of the lapse function (top panel), of the g^x metric compo- 
nent (middle panel) and of the normalized rest-mass density 
(bottom panel). Different lines refer to different times of the 
evolution, with the thick solid line in each panel indicating the 
initial profile and with the thick dashed line corresponding to 
the final timeslice ?At — 0.29 ms; intermediate times (shown 
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FIG. 11. Profiles along the a;-axis of representative metric and 
fluid quantities during the gravitational collapse to a black hole of 
an unstable N — 1, pc — 8.0 x 10""^ relativistic polytrope show- 
ing different snapshots of the time evolution. The top, medium 
and bottom panels show the evolution of the lapse function, of the 
Qxx metric component, and of normalized rest-mass density, respec- 
tively. The thick solid and dashed lines indicate the initial and final 
{t = 0.29 ms) profiles. Intermediate profiles, indicated by thin dot- 
ted ashed lines, are shown every 0.049 ms. 



every 0.049 ms) are indicated with dotted lines. The evolu- 
tion of the lapse function shows the characteristic "collapse 
of the lapse", a distinctive feature of black hole formation. 
The evolution of the g^x metric component and of the rest- 
mass density also clearly exhibit features typical of black hole 
formation, such as the large peak developing in gxx, or the 
continuous increase in the central rest-mass density. 

While the collapse of the lapse is a good indication of the 
formation of a black hole, the formation of an apparent hori- 
zon (the outermost of the trapped surfaces) in the course of the 
simulation is an unambiguous signature of black hole forma- 
tion. An apparent horizon finder based on the fast-flow algo- 
rithm [|T]] was used to detect the appearance of horizons, and 
to calculate the horizon mass. This apparent horizon finder, 
and its validation, is described in Ref. [p^]. 

Fig. |l2| shows the evolution of the horizon mass as a func- 
tion of time. Initially there is no horizon. At a time t ~ 0.21 
ms a black hole forms and an apparent horizon appears. As the 
remaining stellar material continues to accrete onto the newly 
formed black hole its horizon mass increases, finally levelling 
off, until about t = 0.27 ms. The subsequent growth of the 
horizon mass is the result of the increasing error due to grid 
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FIG. 12. Horizon Mass as a function of time. A black hole is 
formed at f = 0.21 ms and the horizon mass then starts to increase 
as a result of accretion. 



Stretching - the radial metric function develops a sharp peak 
which cannot be resolved adequately. 
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FIG. 13. Profiles of the (normalized) rest-mass density along the 
a;-axis and z-axis at two coordinate times, t — (solid lines) and 
t — 3.78 ms (dashed lines), corresponding to three rotational peri- 
ods (P). The star is a A'^ = 1, pc = 1-28 x 10^"^ polytrope rotating 
at 92% of the mass-shedding limit. The simulation has been per- 
formed only in the volume above the {x,y) plane which is covered 
with 129 X 129 x 66 gridpoints. 



V. RAPIDLY ROTATING RELATIVISTIC STARS 

A. Stationary equilibrium models 

The long-term evolution of rapidly rotating, stable equilib- 
rium relativistic stars represents a much more demanding test 
for a numerical code. In this case, in fact, the use of a non-zero 
shift vector is strictly necessary and this, in turn, involves the 
testing of parts of the code that are not involved in the evolu- 
tion of a non-rotating stellar model. The initial data used here 
are numerical solutions describing general relativistic station- 
ary and axisymmetric equilibrium models rotating uniformly 
with angular velocity il. The models are constructed with the 
rns code [ ]53| , ^4| ] (see Ref. [^sj] for a recent review of rotating 
stars in relativity) which provides the four metric potentials i^, 
B, /i, and uj needed to describe the spacetime with line ele- 
ment 



I- r^dO^ 



2v 1 ■ 2 

r sm 



Q{d^ - ujdtf 



(19) 



In the nonrotating limit, the above metric reduces to the metric 
of a static, spherically symmetric spacetime in isotropic coor- 
dinates. A rotating model is uniquely determined upon spec- 
ification of the EOS and two parameters, such as the central 
rest-mass density and the ratio of the polar to the equatorial 
coordinate radii (axes ratio). 

Using the standard Jacobian transformations between the 
spherical polar coordinates (r, 0, (/i) and the Cartesian coor- 



dinates (x, y, z), the initial data for a rotating star are trans- 
formed to Cartesian coordinates. Convergence tests of the 
initial data on the Cartesian grid at various resolutions, show 
that the Hamiltonian and momentum constraints converge at 
second-order everywhere except at the surface of the star, 
where some high-frequency noise is present. This noise is 
due to Gibbs phenomena at the surface of the star, which are 
inherent to the method p^ ] used in the construction of the 
2D initial data (see the relevant discussion in Ref. |]54|]). To 
our knowledge, all currently available methods for construct- 
ing initial data describing rotating relativistic stars suffer from 
some kind of Gibbs phenomena at the surface of the star, 
with the only exception being a recent multi-domain spec- 
tral method that uses surface-adapted coordinates The 
high-frequency noise does not appear to affect the long-term 
evolution of the initial data at the grid resolutions employed 
in our simulations. The evolution is carried out up to several 
rotational periods, using the shift 3-vector obtained from the 
solution of the stationary problem, which we do not evolve in 
time. 

We have evolved models at various rotation rates and for 
several polytropic EOS, all showing similar long-term be- 
haviour and convergence. Hereafter we will focus on a iV = 
1 polytropic model, rotating at 92% of the allowed mass- 
shedding limit for a uniformly rotating star with the same cen- 
tral rest-mass density. In particular, we have chosen a stellar 
model with the same centra l rest-mass density as the nonrotat- 
ing model of Section [V A and which is significantly flattened 
by the rapid rotation (the polar coordinate radius is only 70 % 
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FIG. 14. Profile of the metric component Qxx along the a;-axis 
and z-axis at two different coordinate times, for the same evolution 
shown in Fig. |^. 

of the equatorial coordinate radius). 

Similarly to what is observed in the numerical evolution of 
nonrotating stars, the truncation errors trigger, in a rapidly ro- 
tating star, oscillations that are quasi-radial. As a result, the 
rotating star pulsates mainly in its fundamental quasi-radial 
mode and, during the long-term evolution, its central rest- 
mass density drifts towards higher values. Also in this case, 
both the amplitude of the pulsations and the central density 
growth rate converge to zero at nearly second-order with in- 
creasing grid resolution. 

Our simulations have been performed only in the volume 
above the (x, y) plane which is covered with 129 x 129 x 66 
gridpoints. At such grid resolutions, we have been able to 
evolve a stationary rapidly rotating relativistic star for three 
complete rotational periods, before the numerical solution 
departs noticeably from the initial configuration. Note that 
much longer evolution times (more than an order of magni- 
tude longer and essentially limited by the time available) can 
be achieved if the spacetime is held fixed and only the hy- 
drodynamical equations in a curved background are evolved. 
This has been demonstrated recently in Ref. |]l7|], with a code 
based on the one used in the present paper and in which a 
third-order Piecewise Parabolic Method (PPM) [^] was used 
for the hydrodynamical evolution and applied to the study of 
nonlinear /--modes in rapidly rotating relativistic stars and the 
occuiTence of differential of a kinematical differential rota- 
tion [|9ll (see Ref. [ |60| , ^ for a recent review on the r-mode 
instability). While our cuiTent second-order TVD method 
with the MC limiter is not as accurate (for the same grid reso- 
lution) as the third-order PPM method, it has, nevertheless, a 
very good accuracy, significantly better than that of the min- 
mod limiter. 

Results of our simulations of rapidly-rotating stars are plot- 



ted in Figs. |13|-|15} In particular. Fig. |13| shows the (normal- 
ized) rest-mass density along the x and z axes at two coor- 
dinate times, t = (solid lines) and t — 3.78 ms (dashed 
lines), with the latter coiTesponding to three rotational peri- 
ods. The outer boundary of the grid is placed at about twice 
the equatorial radius. After three rotational periods, the rest- 
mass density profile is still very close to the initial one. Sim- 
ilarly, Fig. |l^ shows the metric component g^x along the x 
and z axes at the same coordinate times of Fig. Again, 
the change in g^x is minimal and only near the stellar surface 
can one observe a noticeable difference (the eiTor there grows 
faster, due to the Gibbs phenomenon in the initial data). 
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FIG. 15. The velocity component along the x-axis at two 
different coordinate times, for the same evolution as in Fig. |lj. 

Besides triggering the appearance of quasi-radial pulsations 
and the secular increase in the central rest-mass density, the 
truncation errors also induce the formation of a local maxi- 
mum at the stellar surface for the evolved "momentum" vari- 
able Sj [cf. Eq. The existence of this local extremum 
reduces, at the surface of the rotating star, the order of our 
TVD schemes to first-order only. As a result, the angular mo- 
mentum profile at the surface gradually drifts away from the 
initial uniformly rotating one, with the rate of convergence of 
this drift being only first-order with increasing grid resolution. 
We emphasize, however, that this is only a local effect: every- 
where else inside the star, the angular momentum evolution is 
second-order accurate. Fig. |l5| shows the velocity component 
yy along the x-axis at the same coordinate times of Fig. ISjand 
p^ . Alternative evolution schemes based on third-order meth- 
ods have been shown to have a smaller truncation error at the 
surface of the star, both for 2D and 3D evolutions of the same 
initial data | p3| , p^ , at least in the Cowling approximation. 

Note that plotting the velocity profile as in Fig. [ij allows 
one to ascertain the accuracy in the preservation of the veloc- 
ity field. Isocontours or vector plots of the velocity field can, 
in fact, easily mask the secular evolution shown in Fig. p3l We 
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also note that the variable evolved in the code is not the ro- 
tational velocity, but a corresponding momentum component 
which depends on the local rest-mass [cf. Eq.( ^]. The er- 
ror in the rotational velocity near the surface is therefore also 
influenced by the small value of the rest-mass density in that 
region. 



B. Quasi-radial modes of rapidly rotating relativistic stars 

The quasi-radial pulsations of rotating neutron stars are a 
potential source of detectable gravitational waves and could 
be excited in various astrophysical scenarios, such as a ro- 
tating core-collapse, a core-quake in a rotating neutron star 
(due to a large phase-transition in the equation of state) or 
the formation of a high-mass neutron star in a binary neu- 
tron star merger An observational detection of such pulsa- 
tions would yield valuable information about the equation of 
state of relativistic stars [^2|]. So far, however, the quasi-radial 
modes of rotating relativistic stars have been studied only un- 
der simplifying assumptions such as in the slow-rotation ap- 
proximation 1 63 ,64]] or in the relativistic Cowling approxima- 



tion [ |4-}j| , |65| ] . The spectrum of quasi-radial pulsations in full 
General Relativity has not been solved to date with perturba- 
tion techniques (see Ref. [|5^] for a recent review of the sub- 
ject). 

In this section we take a step forward in the solution of this 
long standing problem in the physics of relativistics stars and 
obtain the first mode-frequencies of rotating stars in full G en 
eral Relativity and rapid rotation. As done in Section 



IV B for 



the radial pulsation of nonrotating stars, we take advantage 
of the very small numerical viscosity of our code to extract 
physically relevant information from the quasi-radial pertur- 
bations induced by truncation errors. The ability to do so 
demonstrates that our current numerical methods are mature 
enough to obtain answers to new problems in the physics of 
relativistics stars. 

TABLE III. Comparison of small-amplitude quasi-radial pulsa- 
tion frequencies obtained with the present 3D code in fixed space- 
time, with frequencies obtained with an independent 2D code. The 
equilibrium model is a A'' = 1 relativistic polytrope rotating at 92% 
of the mass-shedding limit. 



Mode 


Present 3D code 
(kHz) 


2D code 
(kHz) 


Relative Difference 
(%) 


F 


2.468 


2.456 


0.5 


HI 


4.344 


4.357 


0.3 


H2 


6.250 


6.270 


0.3 



Following the approach outlined in Section fV B| , we have 
first computed the quasi-radial mode frequencies from numer- 
ical evolutions of the GRHydro equations in a fixed space- 
time evolution in order to compare with recent results coming 
from an independent 2D nonlinear evolution code [ 65 1 . Ta- 
ble pl| shows the comparison of between the eigenfrequencies 
computed in the Cowling approximation with the 2D code for 



TABLE IV. Quasi-radial pulsation frequencies for a sequence of 
rotating A'^ = 1 polytropes with rotation rates up to 97% of the 
mass-shedding limit. The frequencies of the fundamental mode F 
and of the first overtone HI are computed from coupled hydrody- 
namical and spacetime evolutions. The ratio of polar to equatorial 

coordinate radii of the rotating models is also shown. 





Q./Q.K 


F (kHz) 


HI (kHz) 


1.000 


0.000 


1.450 


3.958 


0.950 


0.407 


1.411 


3.852 


0.850 


0.692 


1.350 


3.867 


0.825 


0.789 


1.329 


3.894 


0.775 


0.830 


1.287 


3.953 


0.750 


0.867 


1.265 


4.031 


0.725 


0.899 


1.245 


3.974 


0.700 


0.929 


1.247 


3.887 


0.675 


0.953 


1.209 


3.874 


0.650 


0.974 


1.195 


3.717 



the equilibrium model of the previous Section. Note that the 
newly obtained frequencies differ by less than 0.5%, verifying 
that our code can accurately reproduce them. 

Next, we have computed the quasi-radial frequencies in 
coupled hydrodynamical and spacetime evolutions for rapidly 
rotating stars. As mentioned before, this is a novel study and 
the results obtained cannot be compared with coiTesponding 
results in the literature. To study this, we have carried out 
two types of analysis. Firstly, we have followed the same 
procedure used in the case of a nonrotating star case and ob- 
tained the normalized frequency spectrum of oscillations in- 
duced by the truncation errors. Secondly, we have computed 
the frequency spectrum of oscillations triggered by a small 
but specified perturbation. More precisely, we have intro- 
duced the sa me rad ial perturbation in the rest-mass density 
used in Sect. IV D| to induce collapse: i.e. Acos(7rr/2rp), 
where A = 0.02, r is coordinate distance from the center. 



and is the radial coordinate position of the poles. When 
compared, the results of the two treatments indicate that the 
fundamental mode frequency agrees to within 2%, while the 
HI mode near the mass-shedding limit is probably accurate to 
several percent only (at this resolution). 

To study quasi-radial modes of rapidly rotating relativistic 
stars we have built a sequence of models having the same grid 
resolution, the same equation of state and central rest-mass 
density used in the previous section, varying only the rota- 
tion rate il. The sequence starts with a nonrotating star and 
terminates with a star at 97% of the maximum allowed rota- 
tional frequency Q,k = 0.5363 x 10* s^^ for uniformly rotat- 
ing stars (mass-shedding limit). The results of these simula- 
tions are reported in Table ^ and shown in Fig. |l6[ where the 
frequencies of the lowest two quasi-radial modes are shown. 
Interestingly, the fundamental mode-frequencies (solid lines) 
and their first overtones (dashed lines) show a dependence on 
the increased rotation which is similar to the one observed 
for the corresponding frequencies in the Cowling approxima- 
tion [|6|]. 

In particular, the F-mode frequency decreases monotoni- 
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FIG. 16. Quasi-radial pulsation frequencies for a sequence of ro- 
tating A'^ = 1 polytropes and a number of different rotation rates. 
The frequencies of the fundamental mode F (filled squares) and of 
the first overtone HI (filled circles) are computed from coupled hy- 
drodynamical and spacetime evolutions (solid lines). The sequences 
are also compared with the corresponding results obtained from com- 
putations in the relativistic Cowling approximation. 



cally as the maximum rotation rate is approached. Near the 
mass-shedding limit, the frequency is 18% smaller than the 
frequency of the nonrotating star. The difference between the 
F-mode frequency computed here and the corresponding re- 
sult in the Cowling approximation is nearly constant. Thus, 
one can construct an approximate empirical relation for the 
fundamental quasi-radial frequency of rapidly rotating stars, 
using only the coiTesponding frequency in the Cowling ap- 
proximation, i^cowiing and the frequency of the fundamental 
radial mode in the nonrotating limit, Fn=Q- For the particular 
sequence shown above, the empirical relation reads 



Cowline 



1.246) kHz , 



(20) 



and yields the correct frequencies with an accuracy of better 
than 2% for the most rapidly rotating model. More gener- 
ally, if -Fcowiing.si=o is the frequency of the fundamental ra- 
dial mode in the Cowling approximation, then the empirical 
relation can be written as 



F — Fn=o 



Cowlinf 



-F, 



Cowling, SI— 



(21) 



Such an empirical relation is very useful, as it allows one to 
obtain a good estimate of the fundamental quasi-radial mode 
frequency of rapidly rotating stars by solving the hydrody- 
namical problem in a fixed spacetime, rather than solving the 
much more expensive evolution problem in which the space- 
time and the hydrodynamics are coupled. 

The frequency of the HI mode shows a non-monotonic 
decrease as the mass-shedding limit is approached, depart- 



ing from the behavior obtained in the Cowling approxima- 
tion. The oscillations in the frequency at larger rotation 
rates could be due to "avoided crossings" with frequencies 
of other modes of oscillation (We recall that is refeiTed to as 
"avoided crossing" the typical behaviour shown by two eigen- 
frequency curves which approach smoothly but then depart 
from each other without crossing. At the point of closest ap- 
proach, the properties of the modes on each sequence are ex- 
changed [|66||.). Similar avoided crossings have been observed 
also in the Cowling approximation for higher overtones and 
near the mass-shedding limit (see Refs. | p8| , p5| |). Our results 
indicate therefore that the avoided crossings in a sequence 
of relativistic rotating stars occur for smaller rotation rates 
than predicted by the Cowling approximation. This increases 
the importance of avoided crossings and makes the frequency 
spectrum in rapidly rotating stars more complex than previ- 
ously thought. 



VI. GRAVITATIONAL WAVES FROM A PULSATING STAR 

The ability to extract gravitational wave information from 
simulations of relativistic compact objects is an important 
feature of any 3D General Relativistic hydrodynamics code. 
To assess the ability of our code to extract self-consistent 
and accurate gravitational waveforms we have excited simple 
quadrupolar perturbations in our standard spherical = 1 
polytrope. In particular, on the basis of the angular behavior 
of the £ — 2, /-mode in linear perturbation theory, we have 
introduced in the initial model a perturbation in the velocity 
of the form 



ue{t = 0) = A sin {irr/rs) sin ( 



cos ( 



(22) 



where A — 0.02 is the amplitude of the perturbation and is 
the coordinate radius of the star. 

Following York [p5|], we have then constructed the initial 
data for the perturbed model by solving the constraint equa- 
tions for the unperturbed model with added perturbations and 
then proceeded to evolve this solution in time. As a response 
to the initial perturbations, the star has started a series of 
periodic oscillations, mainly in the fundamental quadrupolar 
mode of oscillation. Other, higher-order modes are also ex- 
cited (and observed) but these are several orders of magnitude 
smaller and play no dynamical role. 

As a consequence of the time-varying mass quadrupolar 
triggered by the oscillations, the perturbed star emits gravita- 
tional waves, which are extracted through a perturbative tech- 
nique discussed in detail in Refs. |]67[-p9|], and in which the 
Zerilli function is expanded in terms of spherical harmonics 
with each component being the solution of an ordinary differ- 
ential equation. 

We plot, in Fig. Ill, the ^ = 2. m = component of the 
Zerilli function ^p'^. The upper panel, in particular, shows 
the waverforms as extracted at — 17.7 km (dotted line) 
and at = 23.6 km (solid line) respectively, with the first 
having been rescaled as r^'^/^ to allow a comparison. The 
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FIG. 17. Gravitation-wave extraction {tp^^'^^) of a perturbed 
nonrotating relativistic star, pulsating mainly in the fundamental 
quadrupolar mode. The top panel shows the reseated waverforms 
as extracted at — 17.7 km (dotted line) and at — 23.6 km 
(solid line). The lower panel, on the other hand, shows the amplitude 
of ■!/'^° extracted at — 23.6 km for two different resolutions of 
64^ (dashed line) and 96"^ gridpoints (solid line), respectively. 

very good agreement between the two waveforms is an indi- 
cation that the gravitational waves have reached their asymp- 
totic waveform. The lower panel, on the other hand, shows the 
amplitude of V'^" extracted at r^, — 23.6 km for two differ- 
ent resolutions of 64^ (dashed line) and 96^ gridpoints (solid 
line), respectively. Note that the gravitational wave signal 
converges to a constant amplitude, as the true gravitational- 
wave damping timescale for this mode is several orders of 
magnitude larger than the total evolution time shown. The 
small decrease in the amplitude observed in our numerical 
evolution is thus entirely due to the effective numerical vis- 
cosity of our scheme. At a resolution of 96^ gridpoints, the 
effective numerical viscosity is sufficiently low to allow for a 
quantitative study of gravitational waves from pulsating stars 
over a timescale of many dynamical times (the largest rela- 
tive numerical error estimated on the basis of the simulations 
presented in Fig. |l^ is 3.3%). 

As done in the previous sections, we have compared the fre- 
quencies derived from our numerical simulations with those 
obtained from perturbative techniques for /-mode oscillations 
of = 1 polytropes. Again in this case the comparison has 
revealed a very good agreement between the two approaches. 
As one would expect, the dominant frequency of the gravita- 
tional waves we extract (1.56 kHz) agrees with the fundamen- 
tal quadrupole /-mode frequency of the star (1.58 kHz) [ p7| ] 
to within 1.3 % (at a resolution of 96^ grid-points). 



VII. CONCLUSIONS 

We have presented results obtained with a 3D general rela- 
tivistic code GR_Astro in a comprehensive study of the long- 
term dynamics of relativistic stars. The code has been built by 
the Washington University /Albert Einstein Institute collabora- 
tion for the NASA Neutron Star Grand Challenge Project [|l2|] 
and is based on the Cactus Computational Toolkit [pj|]. The 
simulations reported here have benefited from several new nu- 
merical strategies that have been implemented in the code and 
that concern both the evolution of the field equations and the 
solution of the hydrodynamical equations. In addition to the 
features of the code discussed in paper I, the present version of 
the code can construct various type of initial data representing 
spherical and rapidly rotating relativistic stars, extract gravi- 
tational waves produced during the simulations and track the 
presence of an apparent horizon when formed. 

All of these improvements have allowed tests and perfor- 
mances well superior to those reported in the companion pa- 
per I. With this improved setup we have shown that our code 
is able to succesfuUy pass stringent long-term evolution tests, 
such as the evolution of both, static and rapidly rotating, sta- 
tionary configurations. We have also considered the evolu- 
tion of relativistic stars unstable to either gravitational col- 
lapse or expansion. In particular, we have shown that unsta- 
ble relativistic stars can, in the course of a numerical evolu- 
tion, expand and migrate to the stable branch of equilibrium 
configurations. As an application of this property, we have 
studied the large-amplitude nonlinear pulsations produced by 
the migration. Nonlinear oscillations are expected to accom- 
pany the formation of a proto-neutron star after a supernova 
core-collapse or after an accretion-induced collapse of a white 
dwarf. 

Particularly significant for their astrophysical application, 
we have investigated the pulsations of both rapidly rotating 
and nonrotating relativistic stars and compared the computed 
frequencies of radial, quasi-radial and quadrupolar oscilla- 
tions with the frequencies obtained from perturbative methods 
or from axisymmetric nonlinear evolutions. We have shown 
that our code reproduces these results with excellent accu- 
racy. As a particularly relevant result, we have obtained the 
first mode-frequencies of rotating stars in full general relativ- 
ity and rapid rotation. A long standing problem, such frequen- 
cies had not been obtained so far by other methods. 

In our view the results discussed in this paper have a dou- 
ble significance. Firstly, they establish the accuracy and reli- 
ability of the numerical techniques employed in our code and 
which, to the best of our knowledge, represent the most accu- 
rate long-term 3D evolutions of relativistic stars available to 
date. Secondly, they show that our current numerical methods 
are mature enough for obtaining answers to new problems in 
the physics of relativistic stars. 
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